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We examine the nature of the transition to the antiferromagnetically ordered state in the half-filled three- 
dimensional Hubbard model using the dual-fermion multiscale approach. Consistent with analytics, in the weak- 
coupling regime we find that spin-flip excitations across the Fermi surface are important, and that the strong 
coupling regime is described by Heisenberg physics. In the intermediate interaction, strong correlation regime 
we find aspects of both local and non-local correlations. We analyze the critical exponents of the transition in 
the strong coupling regime and find them to be consistent with Heisenberg physics down to an interaction of 
U/t = \0. 
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I. INTRODUCTION 

The formation of complex states of matter in quantum 
systems is one of the central subjects of condensed matter 
physics. A particular interest is attached to strongly correlated 
fermionic systems, where an interplay between multiple phys¬ 
ical phenomena leads to a rich set of phases and transitions 
between them, including Mott metal-insulator transitions, su¬ 
perconductivity, and magnetisirPEl. Understanding the micro¬ 
scopic mechanism that triggers a phase transition in these sys¬ 
tems is a task of formidable complexitjEl and generally only 
approximate solutions of basic models can be obtained. 

The prototypical model for electronic correlations is the 
fermionic Hubbard model^. It is a benchmark model for a 
variety of theoretical method^ and cold atom experiment 
In three dimensions and for the particle-hole symmetric case, 
the model has a finite temperature transition to the antiferro- 
rnagnetically ordered state at all values of interaction strength 
The character of this transiti on var ies substantially in 
different parts of the phase diagranPE^l and remains under 
debate. 

At weak interactions, unlike in the lower dimensional 
casJiSl^ the appearance of a magnetic order is attributed to 
thermal spin-flip excitations of electrons across the Eermi 
surfacjn At strong electron repulsion, the emergence of 
magnetic or der is attributed to a spin-spin interaction of lo¬ 
cal momentJ2ini. Away from these limits, analytical results 
are not available and these mechanisms may be coexisting, 
competing, or replaced by different physics. 

Describing the properties in this regime has been attempted 
by employing a set of analytic partial summation meth¬ 
ods, including the unrestricted Hartree Loc k app roximatiorPSI 
the two-particle self- consistent approaclP^EI! stroi^ cou¬ 
pling expansionJSES the spin-fluctuation approadPH, and 
the quantum rotor approacfP^, which partially resolve the 
phase diagram. Numerical approaches to the problem in¬ 
clude the c ouple d cluster expansiorP^and lattice Monte-Carlo 
simulationJiSIlSilll at system sizes up to 10^ sites. 

Simu latio ns using the dynamical mean-field theory 
(DMETj32lII] interpolate between the small- and large- 


interaction regimes and describe the fo rmatio n of local mo¬ 
ments upon the increase in mteractiorP^IlII. The method 
includes all local correlationJ^HIll and captures the Eermi- 
surface nesting mechanism of the transitiorP^lbut does not de¬ 
scribe nonlocal corr elatio ns and predicts phase transitions of 
mean-field characteJ^lIlI! Cluster extensions of DMET^SSESI 
provide goo d estimates for the critical temperaturd^. Green’s 
functions^ni and thermodynamics away from the transitiorPI 
but spatial fluctuations are short-ranged and criticality is not 
accessibld^. 

Recent development of diagrammatic multiscale methods 
allowed the extension of these results by incorporating nonlo¬ 
cal correlation^^SEH A treatment using the dynamical vertex 
approximation (DEA), assuming an anomalous dimension ex¬ 
ponent rj =0, showed critical exponents close to the Heisen¬ 
berg valueJ^ll and described the separability of static and dy¬ 
namic correlations in the spectral properties of the modeffl. 

In this paper we study the antiferromagnetic transition in 
the half-filled, three-dimensional (3D) Hubbard model over a 
wide range of interaction strengths U by means of a diff erent 
multiscale scheme, the dual-fermion (DE) approacfP2M21 xhis 
scheme incorporates DMET as a correlated initial point of 
a perturbation expansion, therefore taking the Eermi-surface 
nesting physics into account. The antiferromagnetic nonlocal 
fluctuations are included through a diagrammatic resumma¬ 
tion of a particle-hole ladder series. The method is referred 
to as the ladder DE approacfl^. It has previously been shown 
to correctly capture criticality in the similar Ealicov-Kimball 
model in different dimensionand together with DFA indi¬ 
cated an exponential reduction in Tr acco rding to the Mermin- 
Wagner theorem in two dimensionsSISl 

We study spin correlation functions, critical exponents, and 
spectral properties of the model. By comparison to DMET re¬ 
sults we describe the impact of nonlocal correlations and spin- 
exchange processes. We show that the magnetic properties 
of the model are characterized by the Eermi-surface nesting 
physics at t//f < 4, by a combination of two distinct mecha¬ 
nisms in the strongly correlated regime 4<t//f<10, and by 
local moment Heisenberg physics at large U/t > 10, where t 
is the hopping amplitude. 
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II. MODEL AND METHOD 


We consider the particle-hole symmetric Hubbard model on 
a 3D simple cubic lattice: 


H = Y^ Ek^k.a + E ^ ("'.t - ^)(«a -\)- 

k.o- 1 ^ ^ 


( 1 ) 


Here lattice sites and lattice momenta are labeled i and 
k respectively, CJ denotes the spin projection, gk = 

— 2 f (cos^;c-|-cosA:,. -l-cos^j.) is the electronic dispersion at mo¬ 
mentum k, and U is the Coulomb repulsion between two 
fermions residing at the same site. We take the hopping f = 1 
as the unit of energy. 

The DF approach was developed in Ref. 03 Here we 
briefly introduce the underlying idea. A detailed derivation 
of the formalism can be found in Ref. [53] The DF approach is 
suited for the problem at hand, because it allows us to address 
the Fermi surface nesting physics, local-moment formation, 
and criticality on the same footing. To achieve this, one starts 
from the DMFT solution of the model, which places a quan¬ 
tum impurity at each lattice site. New, so-called DF degrees of 
freedom are then introduced through a Hubbard-Stratonovich 
transformation that couples lattice and DFs locally. The latter 
interact via the vertex functions of the impurity model and me¬ 
diate the coupling between the impurities, which is absent in 
DMFT. After integrating out the lattice fermions, the resulting 
dual action is 


5 =- E 

(B.k.a 


=^( 0 ) 

'^m.k.cT 


+ ( 2 ) 


where / labels the new fermionic degrees of freedom and or 
labels Matsubara frequencies. The bare propagator ^ = 
~ Sa.a represents the momentum-dependent correc¬ 
tion to the DMFT Green’s function. By virtue of the self- 
consistency condition, — i^LkC^S'k^’ local part of 
the DMFT Green’s function is identical to the Green’s func¬ 
tion of the impurity problem. Here N denotes the volume 
of the system. The interaction U is composed of all re¬ 
ducible, fully antisymmetric n-particle vertex functions of the 

impurity, U[f* Ji] « -\y^^afiSyflafi.?flsky+- ■ - where the 
lowest-order term contains the two-particle vertex function, 
defined as 


ymA=S\ V 3 ' (Cl4c34)imp-.?l.?3(5l2534- 514532 ) 


^2‘^4‘- 


( 3 ) 

Combined indices 1 = {wi, ffi} are used in Eq. (|^ to shorten 
notation. 

The complete solution of problem Q is equivalent to the 
lattice problem, and therefore not tractable. The DF ap¬ 
proach proceeds by constructing a low-order approximation 
in terms of the new variables, which corresponds to a sum¬ 
mation of classes of diagrams in the original variables. The 
noninteracting (Gaussian) ensemble of DFs is equivalent to 
the full DMFT solution of the probletrH Because we are in¬ 
terested in the description of the antiferromagnetic phase tran¬ 
sition we restrict the expansion to ladder-type diagrams in the 


particle-hole channel. The dominance of these diagrams can 
be justified by a power counting argument in 1 /M In two di¬ 
mensions, this set of diagrams introduces critical long-range 
fluctuations destroying the spurious order introduced through 
the mean held and leads to the exponential decay of the spin- 
spin coiTelation^. We note that this “second” approximation 
leads to an overestimation of the critical temperature for high 
interaction strengths but has been found to yield the correct 
critical behaviof^. Other approximations are, in principle, 
possible; for instance, providing feedback to the DMFT hy¬ 
bridization function might lead to an improved estimate of 
the critical temperatureSJ and a consideration of cotTections 
from toher-order vertex functions may yield more precise 
result^. Here, however, we keep the DMFT hybridization 
function unchanged, allowing us to study the impact of non¬ 
local correlations in a controllable way. 

The ladder DF makes the assumption that three-particle 
and higher order terms in the interaction U have a small ef¬ 
fect. Numerical evidence of the validity of this assumption 
has been provided in Ref. 14^ The full fermionic frequency 
dependence of the impurity vertex 7 is retained, while restrict¬ 
ing it to a single bosonic frequency, £2 = 0. This approxima¬ 
tion is similar to the construction of an effective functional for 
paramagnons in the same modeiS!, which has been used to 
describe critical properties of quantum many-body systemJ^. 
We have checked that inclusion of more bosonic frequencies 
in the vertex 7 only leads to minor numerical corrections to 
the transition temperature, and does not affect critical proper¬ 
ties (see Appendix [B|. This supports the observation in Ref. 
l48lon the static nature of nonlocal correlations in the model. 

The DMFT impurity problem is solved numerically using 
the continuous-time quantum Monte Carlo hybridization ex¬ 
pansion methocffS®!. The impurity vertex function is ob¬ 
tained on 160 X 160 fermionic frequencies. We sample the 
Brillouin zone on a grid of N=16 x 16 x 16 points and eval¬ 
uate the local and k-dependent Green’s functions ga and 
Gffl k and the static spin susceptibility in reciprocal space 
= 0 ) and real space (5-5y)(£2 = 0 ), from 
which the critical temperature and exponents are extracted. 
The critical slowing-down close to the transition temperature 
Tc is overcome using an annealing scheme with a gradual de¬ 
crease in temperature. Finite-size effects are eliminated by re¬ 
quiring that the correlation length ^ does not reach L/ 6 , where 
L= 16 is the linear system size in units of the lattice constant. 
A detailed description of the calculation procedure is provided 
in Appendix [a] 


III. RESULTS 


Figure [T] shows the phase boundary of the 3D Hubbard 
model obtained by a range of numerical methods. In the 
weak-coupling regime Tc increases as a function of G; the 
slope is accurately described by the random-phase approx¬ 
imation (RPA; unrestricted Hartree-Fock j^ m odified to ac¬ 
count for quantum fluctuations up to G < In the strong¬ 

coupling regi me a n estimate is given by the Heisenberg limit 
Tc = 3.83/G,IS33!, illustrated by the lower dotted line in Fig. 
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Figure 1. (Color online)Transition temperature Tc as a function of U 
(both in units of r = 1) of the Hubbard model in three dimensions at 
particle-hole symmetry, as obtained by DF (filled dark circles) and 
DMFT [dash-dotted (blue) line]. The comparison curves are ob¬ 
tained by the random phase approximation (RPA) from Refs. IT^ 
ancP^lwith a factor of 1/3 in Tc to account for quantum fluctuations^ 
(dashed black line), dynamical cluster approximation [DCA; filled 
(red) hexaaonsjSl lattice quantum Monte Carlo [QMC; filled (green) 
squares]2£l determinantal diagrammatic Monte Carlo [DDMC; filled 
(blue) trianglespHand dynamical vertex approximation [DFA; filled 
(yellow) diamonds jSl. See also Ref. W\ The combination of RPA 
values at small U with the QMdSl the DC/JSl 

and large-C asymp¬ 
totics is plotted as the thick solid black line. 


[T] In the intermediate-coupling regime, estimates of Tc are 
provided by the lattice quantum Monte CarlcP^ and dynami¬ 
cal cluster approximation (DCAj^, extrapolated to the ther¬ 
modynamic limit, which agree within their respective error 
bars and match the rescaled RPA value at t/ = 4. Results of 
the determinantal diagrammatic Monte Carlo methocP^ agree 
at t/ = 6 and 8 and provide a lower Tc at U <5. The results 
from the ladder DFA methoci^ are close, but nonlocal cor¬ 
rections seem slightly overestimated in the intermediate- to 
weak-coupling regime. In order to highlight the overall trend, 
a combined curve interpolated from RPA values at small U, 
averaged lattice quantum Monte CarlcPSIand DCAP^data, and 
large-t/ asymptotics is shown in black. 

The DMFT results, shown as dashed lines in Fig[T] match 
the RPA values up to t/ « 4. Spatial correlations become im¬ 
portant at t/ >4, where DMFT overestimates the value of Tc. 
At large U DMFT approaches a limit of 6/U, coinciding with 
the mean-field estimatd^ (upper dotted line). 

This paper presents ladder DF results, shown in Fig. as 
dark solid circles. The difference between DF and DMFT re¬ 
sults characterizes the contribution from nonlocal magnetic 
fluctuations. At t/ < 4 the DF method shows no deviation 
from the RPA/DMFT results, confirming the validity of the 
local approach. In the intermediate-coupling regime 4 < t/ < 
10, nonlocal corrections are substantial and DF results for 
the transition temperature are consistent with lattice quantum 
Monte Carlo and DCA within error bars. As U is increased 
further, an improvement on DMFT but deviation of Tc to the 
best estimate and the high-temperature series expansion val¬ 
ues becomes apparent. This is evidence of the relevance of 


diagrams that are not included. 

We mark two important values of U in Fig. [T] The point 
t/s = 4 marks the upper limit of validity of DMFT and hence 
the point where nonlocal correlations become important. Co¬ 
incidentally this point is characterized by a change in concav¬ 
ity of the Tc{U) curve, which does not occur in RPA calcula¬ 
tions. We also add a point Uh — 10, at which the DF Tc begins 
to deviate from the best available estimate. This is explained 
in more detail in the following. 
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Figure 2. (Color online) Estimate for the density of states at the 
Fermi level, — j3G(jS /2), as a function of the temperature T and inter¬ 
action strength U. Inset: The cut at constant temperature T = 0.355. 

Figure shows an estimate of the density of states at the 
Fermi level A(a) = 0) ~ —j3G(j3/2). The solid line separates 
the isotropic phase from the symmetry broken antiferromag¬ 
netic region. A horizontal cut for T = 0.355 (dashed line) is 
shown in the inset, illustrating a metal-to-insulator crossover. 
For small U <11$ the isotropic phase exhibits Fermi liquid 
metallic behavior. As U is increased the density of states de¬ 
creases and shows pseudo-gap like behavioi® and at t/ ~ 10 
it becomes completely suppressed. We define Un as the point 
where A(0) —/ 0 and a charge gap opens. 





X 


Figure 3. (Color online) (a) Local, nearest neighbor and next-nearest 
neighbor spin-spin correlation function as a function of interaction U 
atT >Tc. The formation of a well-defined magnetic moment occurs 
at large U. (b) Spatial dependence of the spin susceptibility along the 
real x axis at different values oft/ at T > 7/- 


Figure [^a) shows the magnitude of the local, nearest- 
neighbor, and diagonal spin-spin correlators {S^{i)S^{j)){Q. = 
0 ) in proximity to the transition approached from the high- 
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T-Tc 


Figure 5. (Color line) (a) Inverse magnetic susceptibility and (b) 
inverse correlation length as a function of F — Tc at I7/f = 11, 14, 
and 17, plotted on a logarithmic scale. The thick solid line represents 
the power-law dependence with Heisenberg exponent (a) 1.4 and 

(b) V = 0.7 and the dashed line shows the power-law line with the 
Ising critical exponent (a) 1.24 and (b) v = 0.63. 


Figure 4. (Color online) (a) Inverse static magnetic susceptibility at 
q = {K,n,n) and U/t = 0.25;5.0; 17.0 as obtained within the DF 
(filled symbols, dotted lines) and DMFT (open symbols, solid lines) 
methods and (b) corresponding correlation length. Dashed lines 
show the extrapolated low-temperature behavior, whereas dashed 
horizontal lines in (b) indicate the cutoff value for the inverse cor¬ 
relation length =(>IL. 


temperature side, T > Tc- We observe that nonlocal corre¬ 
lations, represented by the spatial dependence of the spin- 
spin correlation function, are present for all values of U and 
strongly increase for U > Us, consistent with the deviation 
in Tc between DMFT and DF shown in Fig. □ Fig. 0b) 
shows the spatial extent of the spin-spin correlations for dif¬ 
ferent values of U. It illustrates the antiferromagnetic nature 
of the model and shows that spin correlations are small for 
U<Us. 

To characterize the influence of nonlocal correlations, we 
examine the temperature dependence of the susceptibilities 
and correlation lengths in Fig. 0 Plotted is the tempera¬ 
ture dependence of the static inverse magnetic susceptibility 
[Fig. 0a)] and the inverse correlation length [Fig. 0b)] at 
q = {n, n, k) and U = 0.25 (left), U = 5 (center), and U = 11 
(right) as obtained by the ladder DF method and DMFT at 
temperatures where the correlation length ^ < L/6. At small 
U DMFT and DF coincide, as expected from Figs. 0an d0 In 
the intermediate-!/ regime, a difference in the transition tem¬ 
perature (as also depicted in Fig. 0 is visible. At large U, the 
deviation of DF from the DMFT data is substantial and a crit¬ 
ical region with different temperature dependence is visible in 
both plots. 

We now turn to the criticality of the phase transition at large 
U. Figure 0a) shows the temperature dependence of the in¬ 
verse magnetic susceptibility ?Aq= {n, n, n) and in the large¬ 
ly regime (shown here for U/t = 17) as a function of 

T — Tc, plotted on a log-log scale. The power-law dependence 
is resolved in 1.5 decades and can be described by a Heisen- 
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Figure 6 . (Color line) Critical exponents 7 (susceptibility; top panel) 
and V (correlation length; bottom panel) of the Hubbard model as ob¬ 
tained by the DF method at 1/ > Uh [filled (red) diamonds]. A com¬ 
parison with the DGA resulP^ and the DF method for the Falicov- 
Kimball modeP^ is provided. The Heisenberg and Ising exponents 
are shown as dotted and dash-dotted lines, respectively. 


berg exponent of 7 = 1.4 ± 0. 0#^. It is substantially different 
from the Ising value rising = 1 .24, plotted by the dashed line in 
the same graph. This clearly shows the malitative difference 
from the mean-field criticality of DMFT^^. Figure0b) shows 
the temperature dependence of the inverse correlation length, 
plotted on a log-log scale. While available data for the fit in 
the critical region in this case are limited to a single decade, 
there is no uncertainty related to the determination of U, as the 
latter is obtained from Fig. 0a). This substantially increases 
the quality of the fit, which reveals the correlation length ex¬ 
ponent of the Heisenberg universality class 7 = 0.7 ±0.01, 
which is again clearly distinguished from the Ising exponent 
Using “ 0.63. 

A combined analysis of these results is presented in Fig. 0 
Shown are the values and fitting errors of the susceptibility ex¬ 
ponent 7 (top) and the correlation length exponent v (bottom) 
as a function of t/ at 1/ > 10. At lower values of U the crit¬ 
ical region with a different power-law divergence of the sus- 
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ceptibility and the correlation length, compared to DMFT, is 
not numerically accessible. The inhnite-t/ Heisenberg limit 
is also plotted, as the dashed lind^, the Ising result is the 
dot-dashed line. Results from the ladder DTA (two data sets) 
methocPS are shown by (green) pentagons, DF results for the 
3D Falicov-Kimball modei^, obtained using the same fitting 
procedure, are plotted by hlled (blue) circles, and our results 
for the 3D Hubbard model are depicted by (red) diamonds. 
We hnd that, above an interaction strength of t/a = 10, our 
numerical estimates are consistent with the Heisenberg val¬ 
ues. 


IV. DISCUSSION 

Our results provide a straightforward picture of the mag¬ 
netism of the half-filled 3D Hubbard model. At small val¬ 
ues of U < Us the spin-spin interaction between local mo¬ 
ments is not important, as indicated in Fig. The quantita¬ 
tive agreement between the exact transition temperature and 
the estimate, given by the DMFT and DF in Fig. suggests 
that nonlocal corrections to the DMFT Tc are small and hence 
the phase transition to the antiferromagnetic phase is well de¬ 
scribed by a local theory. 

The increase in U leads to a gradual increase of the rel¬ 
evance of nonlocal fluctuations. This is not captured by 
the Fermi-surface nesting mechanism. In the intermediate¬ 
coupling regime Us '^U < Uh, Tc continues to increase (in¬ 
dicative of a Fermi-surface nesting contribution to Tc) and 
nonlocal spin-spin correlations become non-negligible (in¬ 
dicative of Heisenberg physics), pointing to the coexistence 
of both scenarios, revealed in the decrease in the local den¬ 
sity of states. This shows some similarities to the D —i' 
case, obtained by the DCa!^ and DMFT^^ methods, where 
the hnite-temperature transition to the magnetic state is also 
present at all values of U and is not destroyed by Goldstone 
bosons. At f/ = I/h ~ 10 the density of states at the Fermi 
level vanishes, as shown in Fig. Further increase in U leads 
to a decrease in U in accordance with the Heisenberg picture 
and in contradiction to the Fermi-surface nesting picture. 

Unlike the Tc, the critical exponents are universal and, as 
shown in Fig. are recovered by the DF method. The value 
of the exponents match the Heisenberg values precisely at t/ > 
Uh, which conhrms the validity of the Heisenberg picture in 
this parameter regime. 

At first glance it seems surprising that the DF method pro¬ 
vides a correct set of exponents for different models and dif¬ 
ferent dimensions (see also Ref. [50] for an analysis of the 
Falicov-Kimball model), as, in general, an unbiased resum¬ 
mation of multiple competing channels is required and is nor¬ 
mally done using a renormalization-group procedure. When 
only a single instability exists in the system a renormalization- 
group summation results in a ladder set of diagramJ^U. The 
fact that the correct values for the critical exponents are repro¬ 
duced by the DF method indicates that only a single instability 
exists. This points to the importance of the correla ted D MFT 
starting point, which includes all local correlationJ^Hini Val¬ 
idation of this conjecture by a renormalization-group proce¬ 


dure, formulated in the space of dfP, and comparison with 
more involved approximations beyond ladder diagrams, such 
as parquet diagram 

is a direction for future research. 


V. CONCLUSIONS 

We have studied the antiferromagnetic transition in the 
particle-hole symmetric 3D Hubbard model by the ladder 
DF approach and quantihed the relative importance of lo¬ 
cal and nonlocal correlations. We have characterized differ¬ 
ent regimes of the spin-ordering phase transition and inferred 
the physical mechanisms behind them. At t/ < Us « 4, the 
physics is consistent with nesting of the Fermi surface. At 
Us ^ U < Uh ^ 10, both Fermi-surface nesting and spin¬ 
ordering are important. This is the strongly correlated metal¬ 
lic regime, in which Tc increases as a function of U. Spin 
fluctuations reduce the density of states on the Fermi level 
and the position of the critical temperature. However, the 
transition is still governed by the nesting of the Fermi sur¬ 
face. As U > Uh, a charge gap is present and nonlocal cor¬ 
relations become dominant. They manifest themselves in the 
spin-exchange mechanism of the ordering transition, resulting 
in a decrease in transition temperature upon the increase in U. 
The transition is described by Heisenberg critical exponents. 

The ladder DF approach employed here provides a mul¬ 
tiscale description, accounting for correlations on different 
length scales. It may serve as a tool for studying phase transi¬ 
tions and criticality in strongly correlated systems, where the 
application of other numerical methods, such as the quantum 
Monte Carlo or the DCA, is difficult. 
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APPENDIX A: CALCULATION PROCEDURE 

To obtain the building blocks for the dual-perturbation 
theory we solve the impurity problem numerically using 
the continuous-time quantum Monte Carlo hybridization 
expansiorl^. We sample the single particle Green’s function 
on 256 and the impurity vertex functions on 160 fermionic 
and 1 bosonic Matsubara frequency. 
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The Brillouin zone is sampled on a grid ofA^=16xl6xl6 
points. The set of DF equations originates from a skeleton di¬ 
agram expansion of the dual self-energy Ea k- In the ladder 
DF approach, which is used in this work, all diagrams taken 
into account are of the particle-hole ladder type. Summing up 
all these ladder diagrams is most conveniently done by solv¬ 
ing the Bethe-Salpeter equations for the charge and spin com¬ 
ponents of the renormalized particle-hole vertex defined by 
pch/sp _ ptttt pttii and then plugging this solution into the 
Schwinger-Dyson equation. The Bethe-Salpeter equations for 
the charge and spin components of the renormalized particle- 
hole vertex read 


psp/ch _ .^p/ch I 7-T/p/eh y(0) psp/ch 

^ (Bffl'n.q ‘mca'Q. ' ^ ‘coa)"Q. ^ m"a)'n.q 

CD" 

^in.q = ^ JI^^T^CD+n.k+q, 


(A.l) 

(A.2) 


component of the renormalized vertex F^P via 

_y(0) I V' F*P Y^®^ (A St 

AO.k — ^n.k ]\f ^CDO.k^ I 


Note that this relation holds for lattice fermions and dual 
fermions regardless, as long as the corresponding bubble 
or and vertex F®? or f®P are used. However, due to the 
correspondence between Green’s functions of dual and lattice 
fermions at the single- and two-particle level, respectively (for 
details see Ref.Ei), it is possible to relate the spin susceptibil¬ 
ity X'^P for lattice fermions to th e spi n component of the renor¬ 
malized dual vertex f *P if Eq. (A.5 1 is rewritten as 


= (A. 6 ) 

= ~ p; y^, G(3,kGn)+n,k+af^(D.k^cD+fl.k+a ■ (A.7) 

^ k 


where particle-hole notation is adopted and capital Greek let¬ 
ters a re us ed to label bosonic Matsubara frequencies. Equa¬ 
tion (A.2 1 defines the (dual) particle-hole bubble and 
labels the spin and charge components of the impurity ver¬ 
tex. After solving Eq. (A. 1 1 , either by inversion or iteratively, 
the renormalized vertex together with the Schwinger-Dyson 
equation yields the dual self-energy: 



^ H lffl*’(a'a^CD'+n,k+qAffl'a,q 

m'.n.q 

^ 52 lffl^(a'a^®'+n,k+qA(B'n,q 
CD'.tl.q 


■pch 

^ m'con.q 


pSp 

^ m'o)Q..q 


iyCh 
2 ' CD' o£2,q 

ly^p 

2 ‘(n'mQ.q 

(A.3) 


Together with the Dyson equation for DEs [G^k] ^ = 


=.( 0 ) ■ 
'^CD.k 


— k Eqs. (A.l I and (A.31 form a set of nonlin¬ 


ear equations that can be solved self-consistently until conver¬ 
gence in G(B,k is achieved, ensuring the approximation to be 
conserving in the Baym-Kadanoff sense, since all diagrams 
taken into account are dressed skeletons. The solution (fixed 
point) of the ladder DE equations is obtained by an iterative 
solution scheme up to a certain accuracy, which can be esti¬ 


mated by e = k l^a k*^ “ ^ffl.kl’ where g|^]^ is the result 
for the Green’s function obtained in the nth iteration step (here 
spin indices are omitted). A mixing between iterations as low 
as E = 0.05 is used because of the critical slowing-down at 
the phase transition. Resulting dual Green’s functions GcD.k 
are in the convergence radius of the series summation if the 
eigenvalue problem. 


yz 7(B(B'aGo',kGs,/+n,k+q0cD' = ^<pa>, (A.4) 

k.(B' 


yields a leading eigenvalue Amax < 1 • In this case the underly¬ 
ing Bethe-Salpeter equation (BSE) can be solved by inversion. 
Note that due to the connection between the vertex and the 
susceptibility, a leading eigenvalue X,nax = 1 implies a diver¬ 
gent susceptibility. Eor leading eigenvalues below this thresh¬ 
old the spin susceptibility is readily obtained from the spin 


In the equations above the asterisk is a shorthand notation 
to replace the convolution in Eq. ( |A.5| l and the bubble has 
been replaced by X in the second term only. Hence the self- 
consistent solution of Eq. (A. 1 1 allows one to account for the 
influence of spatial correlations onto single and two-particle 
quantities. 

Achieving convergence close to the phase transition is a 
delicate task, as one is confronted with the problem that the 
initial guess provided by is not within the convergence 
radius of the Bethe-Salpeter equation, while the solution is. 
One way to address this problem is to cut off the eigenvalues 
.at 1 during the iterations until the Green’s function lies within 
the convergence radiuJ^U. Here we instead performed suc¬ 
cessive calculations at decreasing temperatures Tq > ... > Ti, 
where we choose our starting point Tq > 
temperature is high enough for g[^]^ to be a sufficient initial 


guess. After obtaining a converged solution G® k at temper¬ 
ature Tm, we interpolate this solution to the Matsubara fre¬ 
quency grid defined by (0„ = ( 2 n -f 1 )for every k to ob¬ 
tain a suitable initial guess for the iterative solution of the DE 
equations at temperature T,n+\. The temperature is then low¬ 
ered with every step, allowing for much faster convergence 
since the information about the spatial modulation of the self¬ 
energy is retained in the initial guess. 

We perform temperature scans, as described above, at dif¬ 
ferent fixed values of the Hubbard U ranging from 0.25 to 17 
in units of the hopping. The correlation length is obtained by 
fitting Xfii^ = c(l + 2 (kjc — to the inverse spin suscep¬ 
tibility = 0,kx,ky = 7Z,k^ = 7z) along the k_j-axis, 

with fit parameterize and \. We determine the critical expo¬ 
nents 7 and V and the critical temperature by fitting a power 
law, x '^K^) = '^(7 ~ to the inverse spin susceptibility 

X^P (n = 0, k = R)^^ and the corresponding correlation length 
= b{T — Tc)^'' with fixed Tc for ^ respectively. The 
fits are performed within the so-called critical region, which 
is the temperature interval [Tj,, 7/] where the behavior of x and 
^ deviates from the high-temperature mean-field dependence. 
While Ti is given by the temperature at which the correlation 
length ^ reaches 1 /6 of the linear system size L, it is more 
difficult to correctly determine the upper boundary Tu of this 
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critical region, which is crucial to avoid ambiguity in the fit. In 
order to do this we employ the procedure described in Ref. |47] 
and used in Ref. |50]in the text and first determine the mean- 
field asymptotics of the inverse susceptibility by fitting a lin¬ 
ear function, j*™, to the high-temperature regime. The upper 
boundary of the critical region Tu is determined as the temper¬ 
ature at which the ratio r = exceeds a certain value 

re- This way we can identify Tu as the temperature where the 
spin susceptibility exits the linear mean-field regime. To elim¬ 
inate the arbitrariness introduced by a specific choice of rc, Tu 
is determined for a set of values of = 0.07-^0.15. 



Figure Al. (Color online) Top: Inverse magnetic susceptibility at 
({= for U/l = 11,14,17 as a function oiT — Tc- Circles, DF 

values; dashed line, linear high-temperature mean-field fit; bold solid 
line, fit of critical exponent in the critical region. Bottom: Residue 
values of the susceptibility fit plotted versus the temperature. 


2.5 
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Figure A2. (Color online) Inverse correlation length (top) and residue 
values of the power-law fit (bottom) at (7/f = 11,14,17 as a function 
of r — Tc- The critical region, obtained from the data in Fig. |A1[ is 
illustrated by the bold line. 


The final result for the critical exponents 7 , v and tempera¬ 
ture Tc is taken as the mean of the obtained fit parameters. The 
error bars are given by a quadrature sum of fit errors and aver¬ 
aging. As shown in the bottom panels in Figs. Al and A2 the 
quality of the fit is satisfactory, showing small temperature- 
independent residues. 


APPENDIX B: IMPACT OF FREQUENCY DEPENDENCE 
OF THE VERTEX 
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Figure A3. (Color online) Same as Fig. |A1| for values of 

U/t = 11,12,13 comparing the static approximation to the DF re¬ 
sults at more than one bosonic frequency. Open symbols, static- 
approximation ladder DF approximation results; filled symbols, 
LDFA results using 120 bosonic frequencies. 


Figure [B| shows a comparison of the inverse magnetic sus¬ 
ceptibilities obtained with the impurity vertex being sam¬ 
pled over 1 (static) or 120 (nonstatic) bosonic Matsubara fre¬ 
quencies, respectively. Two curves are separated by a con¬ 
stant offset, yielding a slightly lower value for the critical 
temperature in the nonstatic case. However both curves ex¬ 
hibit the same slope. Since we estimate the critical expo¬ 
nent 7 from a power-law fit to the inverse spin susceptibil¬ 
ity (T) = a{T — Tc)^^, the static approximation will yield 
the same exponents as the nonstatic estimate since the con¬ 
stant offset is absorbed in the fit parameter 7),. Therefore the 
choice of the number of bosonic Matsubara frequencies does 
not affect the observed scaling quantities. Our approximation 
only takes low-energetic two-particle excitations into account, 
which seem to provide the leading divergent contributions to 
the static magnetic susceptibility and therefore correctly cap¬ 
ture the criticality at the classical phase transition. In order 
to make a quantitative comparison of this approximation in a 
more general case, such as away from particle-hole symmetry, 
it is necessary to include more bosonic frequencies. 
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